!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine bnd_grpr_2d (ro, mx, my, mz, bx, by, bz, en, gx, gy, x, y, gm, nx, ny, ng, dir, sides)
!==================================================================================================================================|
!  sides = 'l' for lower boundary
!  sides = 'h' for higher boundary
!  sides = 'b' for both
!----------------------------------------------------------------------------------------------------------------------------------|

   implicit none

   character(1), intent(in) ::  dir, sides
   integer(4), intent(in) :: ng, nx, ny
   real(8), intent(in) :: gm
   real(8), dimension(nx), intent(in) :: x
   real(8), dimension(ny), intent(in) :: y
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy

   integer(4) :: i, j, ngp1, ngp2, nxm1, nxm0, nym1, nym0, bnd_base, bnd_next
   real(8) :: pr, v2, b2, gmm1
   
!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0
   ngp1 = 1 + ng
   ngp2 = 2 + ng

   if (dir .eq. 'x') then

      if (sides .eq. 'l' .or. sides .eq. 'b') then
         do j = 1, ny
         do i = 1, ng
            bnd_next = ngp1 - i
            bnd_base = ngp2 - i
            v2 = mx(bnd_base, j) * mx(bnd_base, j) + my(bnd_base, j) * my(bnd_base, j) + mz(bnd_base, j) * mz(bnd_base, j)
            b2 = bx(bnd_base, j) * bx(bnd_base, j) + by(bnd_base, j) * by(bnd_base, j) + bz(bnd_base, j) * bz(bnd_base, j)
	         pr = (en(bnd_base, j) - v2 / ro(bnd_base, j) / 2.0d0 - b2 / 2.0d0) * gmm1
            pr = pr + (ro(bnd_next, j) * gx(bnd_next, j) + ro(bnd_base, j) * gx(bnd_base, j)) / 2.0d0 * (x(bnd_next) - x(bnd_base))
            v2 = mx(bnd_next, j) * mx(bnd_next, j) + my(bnd_next, j) * my(bnd_next, j) + mz(bnd_next, j) * mz(bnd_next, j)
            b2 = bx(bnd_next, j) * bx(bnd_next, j) + by(bnd_next, j) * by(bnd_next, j) + bz(bnd_next, j) * bz(bnd_next, j)
            en(bnd_next, j) = v2 / ro(bnd_next, j) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
         enddo
         enddo
      endif

      nxm0 = nx - ng
      nxm1 = nx - ng - 1

      if (sides .eq. 'h' .or. sides .eq. 'b') then
         do j = 1, ny
         do i = 1, ng
            bnd_next = nxm0 + i
            bnd_base = nxm1 + i
            v2 = mx(bnd_base, j) * mx(bnd_base, j) + my(bnd_base, j) * my(bnd_base, j) + mz(bnd_base, j) * mz(bnd_base, j)
            b2 = bx(bnd_base, j) * bx(bnd_base, j) + by(bnd_base, j) * by(bnd_base, j) + bz(bnd_base, j) * bz(bnd_base, j)
	         pr = (en(bnd_base, j) - v2 / ro(bnd_base, j) / 2.0d0 - b2 / 2.0d0) * gmm1
            pr = pr + (ro(bnd_next, j) * gx(bnd_next, j) + ro(bnd_base, j) * gx(bnd_base, j)) / 2.0d0 * (x(bnd_next) - x(bnd_base))
            v2 = mx(bnd_next, j) * mx(bnd_next, j) + my(bnd_next, j) * my(bnd_next, j) + mz(bnd_next, j) * mz(bnd_next, j)
            b2 = bx(bnd_next, j) * bx(bnd_next, j) + by(bnd_next, j) * by(bnd_next, j) + bz(bnd_next, j) * bz(bnd_next, j)
            en(bnd_next, j) = v2 / ro(bnd_next, j) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
         enddo
         enddo
      endif

   elseif (dir .eq. 'y') then

      if (sides .eq. 'l' .or. sides .eq. 'b') then
         do j = 1, ng
            bnd_next = ngp1 - j
            bnd_base = ngp2 - j
            do i = 1, nx
               v2 = mx(i, bnd_base) * mx(i, bnd_base) + my(i, bnd_base) * my(i, bnd_base) + mz(i, bnd_base) * mz(i, bnd_base)
               b2 = bx(i, bnd_base) * bx(i, bnd_base) + by(i, bnd_base) * by(i, bnd_base) + bz(i, bnd_base) * bz(i, bnd_base)
	            pr = (en(i, bnd_base) - v2 / ro(i, bnd_base) / 2.0d0 - b2 / 2.0d0) * gmm1
               pr = pr + (ro(i, bnd_next) * gy(i, bnd_next) + ro(i, bnd_base) * gy(i, bnd_base)) / 2.0d0 * (y(bnd_next) - y(bnd_base))
               v2 = mx(i, bnd_next) * mx(i, bnd_next) + my(i, bnd_next) * my(i, bnd_next) + mz(i, bnd_next) * mz(i, bnd_next)
               b2 = bx(i, bnd_next) * bx(i, bnd_next) + by(i, bnd_next) * by(i, bnd_next) + bz(i, bnd_next) * bz(i, bnd_next)
               en(i, bnd_next) = v2 / ro(i, bnd_next) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
            enddo
         enddo
      endif

      nym0 = ny - ng
      nym1 = ny - ng - 1

      if (sides .eq. 'h' .or. sides .eq. 'b') then
         do j = 1, ng
            bnd_next = nym0 + j
            bnd_base = nym1 + j
            do i = 1, nx
               v2 = mx(i, bnd_base) * mx(i, bnd_base) + my(i, bnd_base) * my(i, bnd_base) + mz(i, bnd_base) * mz(i, bnd_base)
               b2 = bx(i, bnd_base) * bx(i, bnd_base) + by(i, bnd_base) * by(i, bnd_base) + bz(i, bnd_base) * bz(i, bnd_base)
	            pr = (en(i, bnd_base) - v2 / ro(i, bnd_base) / 2.0d0 - b2 / 2.0d0) * gmm1
               pr = pr + (ro(i, bnd_next) * gy(i, bnd_next) + ro(i, bnd_base) * gy(i, bnd_base)) / 2.0d0 * (y(bnd_next) - y(bnd_base))
               v2 = mx(i, bnd_next) * mx(i, bnd_next) + my(i, bnd_next) * my(i, bnd_next) + mz(i, bnd_next) * mz(i, bnd_next)
               b2 = bx(i, bnd_next) * bx(i, bnd_next) + by(i, bnd_next) * by(i, bnd_next) + bz(i, bnd_next) * bz(i, bnd_next)
               en(i, bnd_next) = v2 / ro(i, bnd_next) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
            enddo
         enddo
      endif

   else

      write (*, *) 'Direction error from bnd_grpr_2d.'
      stop

   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine bnd_grpr_2d
